*****************************************************************************************
* Project: Crisis Management and Attitudes towards Democracy (European Union Politics)  *
* Date: 04 February 2025                                                                *
* Authors:  Lea Heyne & Rubén Ruiz-Rufino                                               *
* Do-file 2: Robustness tests                                                           *
* Data: Eurobaromenter EB93.1                                                           *
*****************************************************************************************

/*
_________________________________________________________________________________________________________________________________________________________
Notes on this dofile
_________________________________________________________________________________________________________________________________________________________
This part of the anlaysis perfroms extensive robustnesss test, mostly following Munoz et al. 2020
and their best practices for unexpected events during survey designs (UESD). We use the EB93.1 to 
conduct a quasi-experiment using the adoption of the Recovery and Resilience Fund on 21 July 2020 
as an unexpected event. 

Note on data: EB93.1 comes from file ZA7649_v2-0-0.dta which can be dowlowaded at:
https://www.gesis.org/en/eurobarometer-data-service/data-and-documentation/standard-special-eb/study-overview/eurobarometer-931-za7649-july-august-2020
_________________________________________________________________________________________________________________________________________________________
*/
clear all

*Macros
global root    "C:\pathway"
global figures "$root\figures\"
global tables  "$root\tables\"
global data    "$root\data\" 

use "$data\eb93.1.dta", clear

*********************
* Ado-files & scheme
*********************
ssc install estout, replace 
ssc install blindschemes, replace 
ssc install addplot, replace
ssc install ebalance, replace
ssc install coefplot, replace
ssc install addplot, replace
net install dm79.pkg, from(http://www.stata.com/stb/stb56)
net install dm49.pkg, from(http://www.stata.com/stb/stb39)
ssc install mmat2tex, replace
ssc install webdoc, replace

set scheme plottigblind

///Data preparation 

*Unify West and East Germany
replace country=4 if country==14

*Country-level variables
g emu=0
replace emu=1 if country==1 | country==2  | country==3 | country==4 | country==5 | country==6 | country==8 | country==11 | country==12 | country==13 | country==16 | country==18
label define EMU 1 "France" 2 "Belgium" 3 "Netherlands" 4 "Germany" 5 "Italy" 6 "Luxembourg" 8 " Ireland" 11 "Greece" 12 "Spain" 13 "Portugal" 16 "Finland" 18 "Austria" 
label value country EMU

keep if emu==1

g bailout=country==8 | country==11 | country==12 | country==13 
lab var bailout "Bailout"

*DV - SWD
rename sd18a swd
replace swd=. if swd==5
label var swd "Satisfaction with democracy"
recode swd (1/2=1 "Satisfied")(3/4=0 "Not satisfied")(else=.),into(swd_bi)


*Clean other relevant variables

gen year=.
replace year=2020 if year==.
rename w11 weight
lab var weight "Euro 15"

rename d11 age
lab var age "Age"
rename d8r1 education
lab var education "Education (years)"
recode education (97/98=.)(11=0)
rename d10 gender
lab var gender "Gender (female)"
rename d15a occupation 
gen paidjob = 1
replace paidjob = 0 if occupation == 1 | occupation == 2 | occupation == 3 | occupation == 4 
label var paidjob "In a paid job"
gen housework = 0
replace housework = 1 if occupation == 1
label var housework "Housework"
gen student = 0
replace student = 1 if occupation == 2
label var student "Student"
gen unemployed = 0
replace unemployed = 1 if occupation == 3
label var unemployed "Unemployed" 
gen retired = 0
replace retired = 1 if occupation == 4
label var retired "Retired"
rename capi_cawi CAPI
label var CAPI "Interview: CAPI"
rename p5 cooperation
label var cooperation "Respondent cooperation (low)"

* Country dummies 
tab country, gen(country_)

label var country_1 "France"
label var country_2 "Belgium"
label var country_3 "Netherlands"
label var country_4 "Germany"
label var country_5 "Italy"
label var country_6 "Luxembourg"
label var country_7 "Ireland"
label var country_8 "Greece"
label var country_9 "Spain"
label var country_10 "Portugal"
label var country_11 "Finland"
label var country_12 "Austria"

* code respondent age 

gen age2 = 0 
replace age2 = 1 if age >= 26

gen age3 = 0 
replace age3 = 1 if age >= 28


*Treatment - Date variable
rename p1 date

*date variable where 1st day of fieldwork is 0 in each country

bys country: egen mindate=min(date)
g datenum=date-mindate

/* Generating the different treatment/control variables used throughout the paper: 

- naive --> All respondents
- twenty --> 20 days before and after the announcement 
- ten --> 10 days before and after the announcement 

Those interviewed on the day of the announcement (July 21) are excluded from all the analyses (see Munoz et al).  
*/ 

* Naive 
gen treatment_naive = . 
replace treatment_naive = 0 if date < 17
replace treatment_naive = 1 if date > 17
tab date treatment_naive, m

* Age reduced
gen treatment_age = . 
replace treatment_age = 0 if date < 17 & age3 == 1
replace treatment_age = 1 if date > 17 & age3 == 1
tab date treatment_age, m

* Twenty
gen treatment_twenty = . 
replace treatment_twenty = 0 if date < 17
replace treatment_twenty = 1 if date > 17 & date < 38
tab date treatment_twenty, m

* ten
gen treatment_ten = . 
replace treatment_ten = 0 if date > 6 & date < 17
replace treatment_ten = 1 if date > 17 & date < 28
tab date treatment_ten, m

* Labelling treatment variables 
label def tr 0"Control" 1"Treatment"
label val treatment_* tr

* Interactions 
g interaction_naive=treatment_naive*bailout
lab var interaction_naive "NextGenEUXBailout"

* Generating time running variable where 0 corresponds to the day after the announcement (July 21)

gen time_zero = date - 17 

replace time_zero = . if date == 17

replace time_zero = time_zero - 1 if time_zero > 0

label var treatment_naive "Treatment group"
label var treatment_twenty "Treatment group"
label var treatment_ten "Treatment group"
label var time_zero "Days"

* Generating alternative time running variable where all days before the announcement are coded as 0 (revised paper)

gen time_zero2 = time_zero
replace time_zero2 = 0 if time_zero < 0
label var time_zero2 "Days after treatment"


*****************************
* Entropy balancing weights
*****************************
ebalance treatment_naive	age education gender paidjob unemployed student retired housework cooperation ///
							,  generate(balance_naive) targets(3)
							
							//why smaller sample? ---> missing values, especially in occupation variables and cooperation
							
* Storing entropy balancing descriptives into matrix to generate table
matrix pre = e(preBal)
matrix post = e(postBal)

* Generating variable that captures the units included in the entropy balancing weights estimation
gen sample_bal = e(sample)

* Applying generated weights and survey set

svyset unique_id [pweight=balance_naive], strata(country)

*****************************************
* Figure 3: Effects of NextGenEU on SWD 
*****************************************

preserve
* Panel A: Naive baseline estimation 
eststo m_1: reg swd_bi i.treatment_naive

eststo m_1b: reg swd_bi i.treatment_naive##c.time_zero
margins, dydx (1.treatment) at(time_zero=(0(1)35)) saving(time1, replace)

eststo m_2: reg swd_bi i.treatment_naive##c.bailout if country!=5 //Italy excluded

local n1 = `e(N)'

* Panel B: Baseline estimation + covariate adjustment + FE
eststo m_3: svy: reg swd_bi i.treatment_naive i.country datenum

eststo m_3b: svy: reg swd_bi i.country i.treatment_naive##c.time_zero
margins, dydx (1.treatment) at(time_zero=(0(1)35)) saving(time2, replace)

eststo m_4: svy: reg swd_bi i.treatment_naive##c.bailout datenum if country!=5 //Italy excluded, no FE

local n3 = `e(N)'

* Panel C: Baseline estimation + covariate adjustment + FE, bailout experience (age adjusted)
eststo m_5: svy: reg swd_bi i.treatment_naive i.country datenum if age3 == 1 

eststo m_5b: svy: reg swd_bi i.country i.treatment_naive##c.time_zero if age3 == 1 
margins, dydx (1.treatment) at(time_zero=(0(1)35)) saving(time3, replace)

eststo m_6: svy: reg swd_bi i.treatment_naive##c.bailout datenum if (age3 == 1 & country!=5) //Italy excluded, no FE

local n5 = `e(N)'

* Generating variable capturing the sample included in the main analyses
gen sample_reg = e(sample)

* Generating Figure 1
local plus = char(177)

coefplot	(m_1, msize(medsmall)) (m_2, msize(medsmall) mcolor(gs5) ciopts(lcolor(gs5 gs5))) || ///
			m_3  m_4  || ///
			m_5  m_6,  ///
			xline(0, lpattern(tight_dot) lwidth(medthick) lcolor(greyscale1)) byopts(row(1)) levels(95 90)	///
			bylabels("(A) Full sample" "(B) Full sample, weights+FE" "(C) Bailout experience, weights+FE") nokey	///
			drop(0.bailout 1.bailout *country datenum) coeflabel(1.treatment_naive = "Treatment group"	///
			1.treatment_naive#c.bailout = "Treatment*Bailout" _cons = "Constant") 	///
			aspect(1.5) scheme(plottigblind)
			
addplot 1: , b1title("Effects on SWD", size(small)) norescaling
addplot 2: , b1title("Effects on SWD") norescaling
addplot 3: , b1title("Effects on SWD") norescaling

graph save "$figures\Experiment1a.gph", replace 
graph export "$figures\Experiment1a.pdf", as(pdf) replace 

* Dropping balance weights
svyset, clear
restore

*Generating Marginal plots

combomarginsplot12 time1, recast(scatter) plot1opts(mcolor(gs3)) ciopt(color(gs3%80)) yline (0, lcolor(greyscale1)) name(time1, replace)
combomarginsplot12 time2, recast(scatter) plot1opts(mcolor(gs3)) ciopt(color(gs3%80)) yline (0, lcolor(greyscale1)) name(time2, replace)
combomarginsplot12 time3, recast(scatter) plot1opts(mcolor(gs3)) ciopt(color(gs3%80)) yline (0, lcolor(greyscale1)) name(time3, replace)

graph combine time1 time2 time3, rows(1)

graph save "$figures\Experiment1b.gph", replace 
graph export "$figures\Experiment1b.pdf", as(pdf) replace 

///Note that the design of this plot was manually adjusted in the stata graph editor to improve the style and presentation after compilation.


*********************
*  ONLINE APPENDIX  *
*********************

******************************************************************************************
* Figure A3: Balance tests: Covariate differences between the treatment and control groups
******************************************************************************************

* Some data management is needed before generating Figure with t-tests for the different control/treatment group comparisons

* Reversing treatment variables so that in means comparisons higher numbers means higher (proportion) in treatment group
foreach var of varlist treatment_*{
	gen r_`var' = `var'
	recode r_`var' (1=0) (0=1)
}


* Rescaling some variables for the display of T-tests
gen age_r = age/10
label var age_r "Age (Years/10)"

gen education_r = education/10 
label var education_r "Education (Years/10)"

* Conducting t-test and storing results in matrix to generate the plot (it is necessary to estimate the confidence intervals "manually") 
foreach tr of varlist r_treatment_*{
	matrix mean = J(1,13,.)
	matrix colnames mean = age_r education_r gender paidjob unemployed student retired housework cooperation 
	matrix CI = J(4,13,.)
	matrix colnames CI = age_r education_r gender paidjob unemployed student retired housework cooperation
	matrix rownames CI = ll95 ul95 ll90 ul90
	local i 0
	foreach var of varlist age_r education_r gender paidjob unemployed student retired housework cooperation{
		quietly: ttest `var', by(`tr') 
		local ++ i 
		local diff = r(mu_1) - r(mu_2) 
		matrix mean[1, `i'] = `diff' 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
		matrix CI[1, `i'] = `ll95' \ `ul95' \ `ll90' \ `ul90'
	}
matrix `tr'_m = mean
matrix `tr'_CI = CI
}


* Generating Figure from results stored in matrices 
local plus = char(177)


coefplot	(matrix(r_treatment_naive_m), xline(0)	///
			ci((r_treatment_naive_CI[1] r_treatment_naive_CI[2]) (r_treatment_naive_CI[3] r_treatment_naive_CI[4]))) ///
			|| (matrix(r_treatment_age_m), xline(0, lpattern(solid))	///
			ci((r_treatment_age_CI[1] r_treatment_age_CI[2]) (r_treatment_age_CI[3] r_treatment_age_CI[4])) ///
			pstyle(p1) if(@ll>-.2&@ul<.2)) ///
			(matrix(r_treatment_age_m), xline(0, lpattern(solid)) ///
			ci((r_treatment_age_CI[1] r_treatment_age_CI[2]) (r_treatment_age_CI[3] r_treatment_age_CI[4])) ///
			pstyle(p1) if(@ll<-.6&@ul>.4) ciopts(recast(pcbarrow))) ///
			, byopts(row(2)) xlabel(-.2(.05).2) xscale(range(-.2 .2)) ///
			transform(* = min(max(@,-.2),.2)) nooffset nokey bylabels("Full sample"  "Bailout experience")

graph save "$figures\Figure_A3.gph", replace 
graph export "$figures\Figure_A3.pdf", as(pdf) replace 

****************************************************************************************************
* Figure A4: Effects of the annoucement for multiple bandwidths: no covariate adjustment
****************************************************************************************************

* Generating matrix that includes the effect estimates and the number of units in the treatment and control groups for multiple bandwidths
local i 0 
matrix results = J(47,8,.)
matrix colnames results = bandwidth coef ll95 ul95 ll90 ul90 n_treat n_contr
forval d = 0/46{
	if `d' == 0{
		local d_neg = -1
	}
	if `d' != 0{
		local d_neg = (`d'+1) * -1 
	}
	gen treatment_band = . 
	replace treatment_band = 1 if time_zero >= 0 & time_zero <= `d'
	replace treatment_band = 0 if time_zero < 0 & time_zero >= `d_neg'
	local bandwidth = `d' + 1 
	quietly: reg swd_bi i.treatment_band
	local ++ i 
	local coef = _b[1.treatment_band]
	local degrees = e(df_r)
	local critical_5 = invttail(`degrees', 0.025)
	local se = _se[1.treatment_band]
	local confvalue_5 = `critical_5' * `se'
	local critical_10 = invttail(`degrees', 0.05)
	local confvalue_10 = `critical_10' * `se'
	local ll95 = `coef' - `confvalue_5'
	local ul95 = `coef' + `confvalue_5'
	local ll90 = `coef' - `confvalue_10'
	local ul90 = `coef' + `confvalue_10'
	qui: sum treatment_band if treatment_band == 1
	local n_treat = r(N)
	qui: sum treatment_band if treatment_band == 0 
	local n_contr = r(N)
	matrix results[`i', 1] = `bandwidth',`coef',`ll95',`ul95',`ll90',`ul90',`n_treat',`n_contr'
	drop treatment_band
} 

* Transforming matrix into variables
svmat results, names(col)

* Calculating total height of histogram bars
gen tot_contr = n_treat + n_contr

* Generating figure
local plus = char(177)

graph twoway	(bar tot_contr bandwidth, fcolor(gs4) ///
				yaxis(2) ylabel(0(2000)12000, labsize(vsmall) axis(2)) ytitle("Number of cases in treatment and control groups", axis(2) size(vsmall)) ///
				yscale(axis(2) alt)) ///
				(rbar tot_contr n_treat bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2) legend(off)) ///
				(function y=0, ra(0 50) lstyle(solid) lcolor(red) lwidth(vthin)) ///
				(rspike ll95 ul95 bandwidth, lwidth(0.2) lcolor(black) yaxis(1) yscale(axis(1) alt)) ///
				(rspike ll90 ul90 bandwidth, lwidth(0.4) lcolor(black) yaxis(1)) ///
				(scatter coef bandwidth, msymbol(O) mcolor(black) yaxis(1) ///
				ytitle("Effects on SWD", size(vsmall)) ylabel(, labsize(vsmall) axis(1)) ///
				xtitle("Bandwidth (+/- days) around event", size(vsmall)) xlabel(0(5)50, labsize(vsmall)) legend(off))

graph save "$figures\Figure_A4.gph", replace
graph export "$figures\Figure_A4.pdf", as(pdf) replace

* Droping variables generated from matrix
drop bandwidth coef ll95 ul95 ll90 ul90 n_treat n_contr tot_contr


**************************************************
* Table A1: Covariate distribution (no weighting)
**************************************************
esttab matrix(pre, fmt(3 3 3 3 3 3)) using "$tables/Table_A1.tex", replace ///
					label mlabels(none) style(tex) ///
					mgroups("Treatment group" "Control group", pattern(1 0 0 1 0 0)) mtitles("Mean" "Variance" "Skewness" "Mean" "Variance" "Skewness")
	
esttab matrix(pre, fmt(3 3 3 3 3 3)) using "$tables/Table_A1.html", replace ///
					label mlabels(none) style(tex) ///
					mgroups("Treatment group" "Control group", pattern(1 0 0 1 0 0)) mtitles("Mean" "Variance" "Skewness" "Mean" "Variance" "Skewness")
								
*************************************************************
* Table A2: Covariate distribution (reweighted control group)
*************************************************************		
esttab matrix(post, fmt(3 3 3 3 3 3)) using "$tables/Table_A2.tex", replace ///
					label mlabels(none) style(tex) ///
					mgroups("Treatment group" "Control group", pattern(1 0 0 1 0 0)) mtitles("Mean" "Variance" "Skewness" "Mean" "Variance" "Skewness")
					
esttab matrix(post, fmt(3 3 3 3 3 3)) using "$tables/Table_A2.html", replace ///
					label mlabels(none) style(tex) ///
					mgroups("Treatment group" "Control group", pattern(1 0 0 1 0 0)) mtitles("Mean" "Variance" "Skewness" "Mean" "Variance" "Skewness")
					

********************************
* Table A3: Placebo treatments.
********************************

* Generating placebo treatment variable at the empirical median of the control group (July 16)
gen placebo = . 

sum date if treatment_twenty == 0 , det

replace placebo = 0 if date <= 11 

replace placebo = 1 if date > 11 & date <= 16

tab placebo

tab date placebo, m

label var placebo "Placebo treatment"
label var date "Fieldwork day"


* Placebo at the empirical median of the control group
eststo p_1: reg swd_bi i.placebo i.country
estadd local ctryFE "Yes"

gen time_zero_p = date - 11
label var time_zero_p days
eststo p_2: reg swd_bi i.placebo##c.time_zero_p i.country
estadd local ctryFE "Yes" 

drop placebo time_zero_p


* Analysis of the relationship between SWD and time during the control period 
gen time_zero_p=time_zero
label var time_zero_p days

eststo p_3: reg swd_bi time_zero_p i.country if treatment_naive == 0 
estadd local ctryFE "Yes"

eststo p_4: reg swd_bi c.time_zero_p##c.time_zero_p i.country if treatment_naive == 0 
estadd local ctryFE "Yes"
 
 
* Generating table
esttab p_* using 	"$tables/Table_A3.tex",  ///
					b(2) se(2) r2   ///
					nobase label replace interaction(" X ")  ///
					nobaselevels  nodiscrete ///
					mlabels(none) drop  (*country*) ///
					scalars("ctryFE Country Fixed-Effects") /// 
					mgroups("Control group"  ///
					"Pre-event period", pattern(1 0 1 0 0 0))
					
esttab p_* using 	"$tables/Table_A3.html",  ///
					b(2) se(2) r2   ///
					nobase label replace interaction(" X ")  ///
					nobaselevels  nodiscrete ///
					mlabels(none) drop  (*country*)  ///
					scalars("ctryFE Country Fixed-Effects") /// 
					mgroups("Control group"  ///
					"Pre-event period", pattern(1 0 1 0 0 0))
					
drop time_zero_p

******************************************************************************
* Figure A5: Average satisfaction with democracy by day in pre-event period
******************************************************************************

* Eliminate outlier (3 interviews on day -16)
drop if time_zero == -16

* Generating variable measuring average satsifaction with government by day 
bysort time_zero: egen avgsatisfaction = mean(swd_bi)

* Generating figure 
graph twoway	(scatter avgsatisfaction time_zero if time_zero < 0,  ysc(r(0 1))  ylabel(0(0.2)1) ) ///
				(lpoly swd_bi time_zero if time_zero < 0) /// 
				(lowess swd_bi time_zero if time_zero < 0) ///
				(lfit swd_bi time_zero if time_zero < 0, lpattern(solid)), /// 
				xtitle("Days until event") ytitle("Satisfaction with democracy") ///
				legend(order(2 "Kernel weighted local polynomial smoothing" 3 "Locally weighted regression smoothing" 4 "Linear Fit") size(vsmall) rows(1) pos(6))
				
graph save "$figures/Figure_A5.gph", replace 
graph export "$figures/Figure_A5.pdf", as(pdf) replace 


********************************************************************
* Figure A6: Falsification tests based on other outcome variables
********************************************************************

recode qa12_3 qa12_1 qa6a_10 qa6a_3 (2=0)(3=.)

recode d70 (2=1)(3/4=0)(5=.)

recode qb8_1 qb8_2 (5=.)

* Re-applying original survey weights 
svyset unique_id [pweight=balance_naive], strata(country)

* Estimating the models: 
	* Trust in legal system 
eststo f_1: svy: reg qa6a_3 i.treatment_naive i.country
eststo f_2: svy: reg qa6a_3 i.treatment_naive##c.bailout if country!=5 //Italy excluded, no FE
local f_1 = `e(N)'

	* Immigration preferences 
eststo f_4: svy: reg qb8_2 i.treatment_naive i.country
eststo f_5: svy: reg qb8_2 i.treatment_naive##c.bailout if country!=5 //Italy excluded, no FE
local f_2 = `e(N)'


* Generating figure
coefplot	(f_1, msize(medsmall)) (f_2, msize(medsmall) mcolor(gs9) ciopts(lcolor(gs9 gs9))) || ///
			f_4 f_5 || , ///
			 xline(0, lpattern(solid)) byopts(row(1)) levels(95 90)	///
			bylabels("(A) Trust Legal System N=7631" "(B) Immigration Preferences N=7456") nokey	///
			coeflabel(1.treatment_naive = "Treatment group" bailout = "Bailout" ///
			1.treatment_naive#c.time_bailout = "Treatment*Bailout"_cons = "Constant") drop (*country) aspect(0.75)


addplot 1: , b1title("Effects on trust legal system", size(small)) norescaling
addplot 2: , b1title("Effects on immigration preferences") norescaling

graph save "$figures/Figure_A6.gph", replace 
graph export "$figures/Figure_A6.pdf", as(pdf) replace 

* Dropping balance weights
svyset, clear

*****************************************************************************
* Figure A2: Google trends for key search terms during the survey fieldwork
*****************************************************************************

use "data\google.dta", clear

graph twoway	(line recoveryfund date, ///
				xline(22117, lpattern(solid)) xlabel(22097 22111 22127 22142 22158, format("%td")) ///
				xtitle("") ytitle("") ti("Recovery Fund")) 
				
graph save rf.gph, replace
				
graph twoway	(line recoveryplan date, ///
				xline(22117, lpattern(solid)) xlabel(22097 22111 22127 22142 22158, format("%td")) ///
				xtitle("") ytitle("") ti("Recovery Plan"))
				
graph save rp.gph, replace

graph twoway	(line nextgen date, ///
				xline(22117, lpattern(solid)) xlabel(22097 22111 22127 22142 22158, format("%td")) ///
				xtitle("") ytitle("") ti("Next Generation EU"))
				
graph save ng.gph, replace

graph combine 	rf.gph ///
				rp.gph ///
				ng.gph, ///
				l1(Relative Google searches in Europe) ///
				b1(Month)
				
				
graph save "$figures/Figure_A2.gph", replace 
graph export "$figures/Figure_A2.pdf", as(pdf) replace 	
